################################################################################
#       Copyright (C) 2010 Michael Yurko <myurko@gmail.com>
#
#This file is part of pyQAP.
#
#pyQAP is free software: you can redistribute it and/or modify
#it under the terms of the GNU General Public License as published by
#the Free Software Foundation, either version 2 of the License, or
#(at your option) any later version.
#
#pyQAP is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with pyQAP.  If not, see <http://www.gnu.org/licenses/>.
################################################################################

import numpy as np
import multiprocessing as mp
import heuristics
'''
Contains helper functions for the heuristics.
'''

def packed_annealing(tuple_of_inputs):
    '''
    Helper function used in the parallel simulated annealing function.
    
    :param tuple_of_inputs: the inputs from parallel annealing
    :type tuple_of_inputs: tuple
    :returns: the result of on annealing
    :rtype: PSolution
    
    Examples:
    
    >>> from pyQAP import *                                                                  
    >>> prob = problems.QAPLIB('nug12')
    >>> packed_annealing((prob, 2, 100))
    Cost: ...
    Partial QAP solution: [...] with 12 placed.
     Placed:12
     Unplaced:[0 0 0 0 0 0 0 0 0 0 0 0]
    >>> prob = problems.QAPLIB('sko100f')
    >>> packed_annealing((prob, 2, 100))
    Cost: ...
    Partial QAP solution: [...] with 100 placed.
     Placed:100
     Unplaced:[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
     0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
    >>> prob = problems.QAPLIB('had16')
    >>> packed_annealing((prob, 2, 100))
    Cost: ...
    Partial QAP solution: [...] with 16 placed.
     Placed:16
     Unplaced:[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
    '''
    problem, mult, steps = tuple_of_inputs
    np.random.seed()
    return heuristics.simulated_annealing(problem, mult,  steps)
    
def parallel_annealing_all(n, threads, problem, mult,  steps):
    '''
    A parallel version of simulated annealing that gives the results of all
    annealings.
    
    Performs simulated annealing in parallel n times with the other 
    given parameters. See simulated_annealing for more info. Returns a list of 
    all the solutions that were found.
    
    :param n: the number of annealings to do
    :type n: int
    :param threads: the number of threads to use
    :type threads: int
    :param problem: the problem to use
    :type problem: Problem
    :param mult: the constant multipluier for the number of iterations to do at
    each temperature
    :type mult: int
    :param steps: the approximate number of temperature steps to preform
    :type steps: int
    
    Examples:
    
    >>> from pyQAP import *
    >>> prob = problems.QAPLIB('nug30')
    >>> vars.verbose = False
    >>> parallel_annealing_all(10, 2, prob, 2, 100)
    [Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0], Partial QAP solution: [...] with 30 placed.
     Placed:30
     Unplaced:[0 ... 0]]

    '''
    pool = mp.Pool(threads)
    psols = pool.map(packed_annealing, ((problem, mult, steps)\
                                         for i in xrange(n)))
    return psols
    
def proj_e(x, a=None, q=None, r=None, b=None):
    '''
    Projection onto the doubly stochastic matrices.
    
    Takes a matrix x and returns the projection onto the doubly stochastic
    matrices. Otherwise, it takes the matrices. Used in gradient_projection.
    
    :param x: the permutation matrix
    :type x: ndarray of doubles
    
    >>> import numpy as np
    >>> np.random.seed(123)                                
    >>> x = np.random.random((4,4))
    >>> e, a, q, r, b = proj_e(x)
    >>> print e                  
    [[ 0.43169687  0.31496729  0.05873023  0.19460562]
     [ 0.19284799  0.19008576  0.55079432  0.06627193]
     [ 0.17003398  0.37481987  0.12893119  0.32621496]
     [ 0.20542116  0.12012708  0.26154426  0.41290749]]
    >>> print a                                        
    [[ 1.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  1.  1.  1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  1.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  1.  1.  1.  1.]
     [ 0.  1.  0.  0.  0.  1.  0.  0.  0.  1.  0.  0.  0.  1.  0.  0.]
     [ 0.  0.  1.  0.  0.  0.  1.  0.  0.  0.  1.  0.  0.  0.  1.  0.]
     [ 0.  0.  0.  1.  0.  0.  0.  1.  0.  0.  0.  1.  0.  0.  0.  1.]]
    >>> np.absolute(q)
    matrix([[  5.00000000e-01,   0.00000000e+00,   0.00000000e+00,
               0.00000000e+00,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01],
            [  5.00000000e-01,   0.00000000e+00,   0.00000000e+00,
               0.00000000e+00,   4.33012702e-01,   1.38777878e-16,
               1.66533454e-16],
            [  5.00000000e-01,   0.00000000e+00,   0.00000000e+00,
               0.00000000e+00,   1.44337567e-01,   4.08248290e-01,
               1.11022302e-16],
            [  5.00000000e-01,   0.00000000e+00,   0.00000000e+00,
               0.00000000e+00,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01],
            [  0.00000000e+00,   5.00000000e-01,   0.00000000e+00,
               0.00000000e+00,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01],
            [  0.00000000e+00,   5.00000000e-01,   0.00000000e+00,
               0.00000000e+00,   4.33012702e-01,   5.55111512e-17,
               5.55111512e-17],
            [  0.00000000e+00,   5.00000000e-01,   0.00000000e+00,
               0.00000000e+00,   1.44337567e-01,   4.08248290e-01,
               2.77555756e-17],
            [  0.00000000e+00,   5.00000000e-01,   0.00000000e+00,
               0.00000000e+00,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01],
            [  0.00000000e+00,   0.00000000e+00,   5.00000000e-01,
               0.00000000e+00,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01],
            [  0.00000000e+00,   0.00000000e+00,   5.00000000e-01,
               0.00000000e+00,   4.33012702e-01,   6.93889390e-17,
               1.11022302e-16],
            [  0.00000000e+00,   0.00000000e+00,   5.00000000e-01,
               0.00000000e+00,   1.44337567e-01,   4.08248290e-01,
               5.55111512e-17],
            [  0.00000000e+00,   0.00000000e+00,   5.00000000e-01,
               0.00000000e+00,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01],
            [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
               5.00000000e-01,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01],
            [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
               5.00000000e-01,   4.33012702e-01,   2.77555756e-17,
               5.55111512e-17],
            [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
               5.00000000e-01,   1.44337567e-01,   4.08248290e-01,
               8.32667268e-17],
            [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
               5.00000000e-01,   1.44337567e-01,   2.04124145e-01,
               3.53553391e-01]])
    >>> np.absolute(r)
    matrix([[ 2.        ,  0.        ,  0.        ,  0.        ,  0.5       ,
              0.5       ,  0.5       ],
            [ 0.        ,  2.        ,  0.        ,  0.        ,  0.5       ,
              0.5       ,  0.5       ],
            [ 0.        ,  0.        ,  2.        ,  0.        ,  0.5       ,
              0.5       ,  0.5       ],
            [ 0.        ,  0.        ,  0.        ,  2.        ,  0.5       ,
              0.5       ,  0.5       ],
            [ 0.        ,  0.        ,  0.        ,  0.        ,  1.73205081,
              0.57735027,  0.57735027],
            [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
              1.63299316,  0.81649658],
            [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
              0.        ,  1.41421356]])
    >>> print b
    [[ 1.]
     [ 1.]
     [ 1.]
     [ 1.]
     [ 1.]
     [ 1.]
     [ 1.]]
    >>> proj_e(e, a, q, r, b)[0]
    matrix([[ 0.43169687,  0.31496729,  0.05873023,  0.19460562],
            [ 0.19284799,  0.19008576,  0.55079432,  0.06627193],
            [ 0.17003398,  0.37481987,  0.12893119,  0.32621496],
            [ 0.20542116,  0.12012708,  0.26154426,  0.41290749]])
    '''
    n = np.size(x,1)
    
    #first time, compute q and r
    if a == None:
        e = np.asmatrix(np.ones((n, 1)))
        a = np.asmatrix(np.zeros((2*n-1, n*n)))
        b = np.asmatrix(np.ones((2*n-1, 1)))
        for i in xrange(n):
            a[i, (i)*n:(i)*n+n] = e.T[:,0]
        for i in xrange(n, 2*n-1):
            for j in xrange(n): 
                a[i, (i-n-3)+(j)*n]=1
        q, r = np.linalg.qr(a.T)
    y=x.copy()
    y.resize(n*n,1)
    #main call after setup
    w = np.linalg.solve(r.T, a*y-b)
    y = y-q*w
    e = y.reshape((n, n))
    return e, a, q, r, b
        
def round_to_perm(x):
    """
    Starting from a matrix close to a permutation matrix, round it to a
    permutation matrix.
    
    :param x: the close to permutation matrix
    :type x: ndarray of doubles
    :returns: a permutation matrix
    :rtype: matrix of zeros
    
    Examples:
    
    >>> import numpy as np
    >>> np.random.seed(123)
    >>> x = np.random.random((4,4))
    >>> x                          
    array([[ 0.69646919,  0.28613933,  0.22685145,  0.55131477],
           [ 0.71946897,  0.42310646,  0.9807642 ,  0.68482974],
           [ 0.4809319 ,  0.39211752,  0.34317802,  0.72904971],
           [ 0.43857224,  0.0596779 ,  0.39804426,  0.73799541]])
    >>> round_to_perm(x)
    matrix([[ 1.,  0.,  0.,  0.],
            [ 0.,  0.,  1.,  0.],
            [ 0.,  1.,  0.,  0.],
            [ 0.,  0.,  0.,  1.]])
    >>> x = np.random.random((8,8))
    >>> x
    array([[ 0.18249173,  0.17545176,  0.53155137,  0.53182759,  0.63440096,
             0.84943179,  0.72445532,  0.61102351],
           [ 0.72244338,  0.32295891,  0.36178866,  0.22826323,  0.29371405,
             0.63097612,  0.09210494,  0.43370117],
           [ 0.43086276,  0.4936851 ,  0.42583029,  0.31226122,  0.42635131,
             0.89338916,  0.94416002,  0.50183668],
           [ 0.62395295,  0.1156184 ,  0.31728548,  0.41482621,  0.86630916,
             0.25045537,  0.48303426,  0.98555979],
           [ 0.51948512,  0.61289453,  0.12062867,  0.8263408 ,  0.60306013,
             0.54506801,  0.34276383,  0.30412079],
           [ 0.41702221,  0.68130077,  0.87545684,  0.51042234,  0.66931378,
             0.58593655,  0.6249035 ,  0.67468905],
           [ 0.84234244,  0.08319499,  0.76368284,  0.24366637,  0.19422296,
             0.57245696,  0.09571252,  0.88532683],
           [ 0.62724897,  0.72341636,  0.01612921,  0.59443188,  0.55678519,
             0.15895964,  0.15307052,  0.69552953]])
    >>> round_to_perm(x)
    matrix([[ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
            [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
            [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.]])
    """
    n = np.size(x, 1)
    x_new = x.copy()
    p = np.asmatrix(np.zeros((n, n)))
    for k in range(n):
        #find max
        max_yet = -float('inf')
        pos = None
        for i in xrange(n):
            for j in xrange(n):
                if max_yet < x_new[i, j]:
                    max_yet = x_new[i, j]
                    pos = (i, j)
        p[pos[0], pos[1]] = 1
        #turn used rows and columns to -inf
        #cols
        for i in xrange(n):
            x_new[pos[0], i] = -float('inf')
        #rows
        for j in xrange(n):
            x_new[j, pos[1]] = -float('inf')
    return p
    
def proj_n(x):
    '''
    Projection onto the non-negative matrices.
    
    If an element is less than zero, then it is set equal to zero.
    
    :param x: a matrix
    :type x: ndarray
    :returns: a non-negative ndarray
    :rtype: ndarray
    
    Examples:
    
    >>> import numpy as np
    >>> np.random.seed(123)
    >>> x = np.random.random((4,4))-.5
    >>> x
    array([[ 0.19646919, -0.21386067, -0.27314855,  0.05131477],
           [ 0.21946897, -0.07689354,  0.4807642 ,  0.18482974],
           [-0.0190681 , -0.10788248, -0.15682198,  0.22904971],
           [-0.06142776, -0.4403221 , -0.10195574,  0.23799541]])
    >>> proj_n(x)
    array([[ 0.19646919,  0.        ,  0.        ,  0.05131477],
           [ 0.21946897,  0.        ,  0.4807642 ,  0.18482974],
           [ 0.        ,  0.        ,  0.        ,  0.22904971],
           [ 0.        ,  0.        ,  0.        ,  0.23799541]])
    >>> x = np.random.random((5,5))-.5
    >>> x
    array([[-0.31750827, -0.32454824,  0.03155137,  0.03182759,  0.13440096],
           [ 0.34943179,  0.22445532,  0.11102351,  0.22244338, -0.17704109],
           [-0.13821134, -0.27173677, -0.20628595,  0.13097612, -0.40789506],
           [-0.06629883, -0.06913724, -0.0063149 , -0.07416971, -0.18773878],
           [-0.07364869,  0.39338916,  0.44416002,  0.00183668,  0.12395295]])
    >>> proj_n(x)
    array([[ 0.        ,  0.        ,  0.03155137,  0.03182759,  0.13440096],
           [ 0.34943179,  0.22445532,  0.11102351,  0.22244338,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.13097612,  0.        ],
           [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ],
           [ 0.        ,  0.39338916,  0.44416002,  0.00183668,  0.12395295]])
    '''
    x[x<0.0]=0
    return x
    
    
def proj_o(x,rev):
    '''
    A projection onto the Orthogonal matrices
    
    Since this is not a convex set, the projection is not unique. Currently uses
    just a QR decomposition. Attempts to coerce a better solution were not
    fullyfimplemented.
    
    :param x: the permutation matrix
    :type x: matrix of doubles
    :param rev: whether the current iteration is odd or not
    :type rev: boolean or int
    
    Examples:
    
    >>> import numpy as np
    >>> np.random.seed(123)
    >>> x = np.random.random((4,4))
    >>> x
    array([[ 0.69646919,  0.28613933,  0.22685145,  0.55131477],
           [ 0.71946897,  0.42310646,  0.9807642 ,  0.68482974],
           [ 0.4809319 ,  0.39211752,  0.34317802,  0.72904971],
           [ 0.43857224,  0.0596779 ,  0.39804426,  0.73799541]])
    >>> np.absolute(proj_o(x,0))
    array([[ 0.58316232,  0.27250708,  0.63968139,  0.42008252],
           [ 0.60242033,  0.25579688,  0.67013808,  0.35010378],
           [ 0.40269027,  0.63148893,  0.29239417,  0.59461579],
           [ 0.367222  ,  0.67935969,  0.23712513,  0.58939807]])
    >>> x = np.random.random((5,5))
    >>> np.absolute(proj_o(x,1))
    array([[ 0.16286111,  0.02631787,  0.91780247,  0.17134089,  0.31790652],
           [ 0.7580585 ,  0.30683639,  0.32886135,  0.20063357,  0.42754546],
           [ 0.32287109,  0.29120927,  0.1780586 ,  0.84755707,  0.24676623],
           [ 0.38704798,  0.03477287,  0.13145143,  0.42234049,  0.80829063],
           [ 0.38048874,  0.9050673 ,  0.02246105,  0.18347749,  0.0437381 ]])
    '''
#    n = np.size(x, 1)
#    if rev:
#        isteps = range(n-1,-1,-1)
#    else:
#        isteps = range(n)
#    q = np.asmatrix(np.zeros((n,n)))
#    for i in isteps:
#        if rev:
#            jsteps = range(n-1,i, -1)
#        else:
#            jsteps = range(0, i)
#        v = x[:,i]
#        for j in jsteps:
#            v=v-v.T*q[:,j]*q[:,j]
#        v=v/np.sqrt(v.T*v)
#        if v[v<0].sum() > v[v>0].sum():
#            v = -v
#        print v
#        q[:, i]=v[:, 0]
    q, r = np.linalg.qr(x)
    return q
    
def rand_perm(n):
    '''
    Returns a random permutation matrix of size n.
    
    The returned matrix is a random permutation matrix of zeros.
    
    :param n: the size of matrix to return
    :type n: int
    
    Examples:
    
    >>> np.random.seed(123)
    >>> rand_perm(5)
    matrix([[ 0.,  1.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  1.,  0.],
            [ 0.,  0.,  0.,  0.,  1.],
            [ 1.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  1.,  0.,  0.]])
    >>> rand_perm(10)                 
    matrix([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
            [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
            [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])
    >>> rand_perm(12)
    matrix([[ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
            [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
            [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
            [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.]])
    >>> rand_perm(8).sum()
    8.0
    >>> rand_perm(50).sum()
    50.0
    '''
    #get random permutation array
    p_array = np.random.permutation(n)
    #create random permutation matrix from p_array
    ret = np.asmatrix(np.zeros((n, n)))
    for i in xrange(n):
        ret[i, p_array[i]] = 1
    return ret
    
def qapobj(x,problem):
    '''
    Return the cost of a given permutation matrix.
    
    :param x: the permutation matrix
    :type x: binary (0-1) matrix
    :param problem: the problem
    :type problem: Problem
    :returns: the cost of the given permutation matrix with the given problem
    :rtype: float
    
    Examples:
    
    >>> from pyQAP import *
    >>> import numpy as np
    >>> np.random.seed(123)
    >>> x = rand_perm(12)  
    >>> prob = problems.QAPLIB('nug12')
    >>> qapobj(x,prob)
    772.0
    >>> x = rand_perm(12)
    >>> qapobj(x,prob)
    918.0
    >>> prob = problems.QAPLIB('had16')
    >>> x = rand_perm(16)
    >>> qapobj(x,prob)
    4278.0
    >>> x = rand_perm(16)
    >>> qapobj(x,prob)
    4280.0
    '''
    return np.trace(problem.f*x*problem.d*x.conj().T)
    
def convert_x_to_array(x):
    '''
    Return the solution array corresponding to the given permutation
    matrix, x.
    
    :param x: the permutation matrix
    :type x: binary (0-1) matrix
    :returns: the corresponding solution array
    :rtype: array
    
    Examples:
    
    >>> from pyQAP import *
    >>> import numpy as np
    >>> np.random.seed(123)
    >>> x = rand_perm(4)
    >>> x
    matrix([[ 0.,  0.,  0.,  1.],
            [ 1.,  0.,  0.,  0.],
            [ 0.,  1.,  0.,  0.],
            [ 0.,  0.,  1.,  0.]])
    >>> convert_x_to_array(x)
    array([3, 0, 1, 2])
    >>> x = rand_perm(10)
    >>> convert_x_to_array(x)
    array([4, 0, 7, 5, 8, 3, 1, 6, 9, 2])
    >>> x = rand_perm(100)
    >>> convert_x_to_array(x)
    array([84, 86, 87, 42, 62, 31, 13, 56, 38,  8,  5,  0, 74, 65,  4, 63, 88,
           28, 41, 70, 81, 33, 29, 57, 53, 23,  9, 54, 24, 19, 95, 35, 90, 72,
           21, 50, 94, 98, 85, 71, 79, 17, 82, 12, 20, 93, 37, 43,  1, 91, 11,
           92, 60, 14, 59, 16, 26, 51, 45,  6, 77, 30, 44, 89, 76, 15, 18, 22,
           10, 58, 75, 64, 69,  3, 40, 34, 27, 52,  7, 48, 61, 99, 66, 39,  2,
           67, 55, 49, 68, 80, 36, 78, 83, 25, 46, 32, 73, 47, 96, 97])
    >>> x = rand_perm(30)
    >>> convert_x_to_array(x)
    array([ 5, 10, 22,  2, 11, 16,  8, 23,  9,  0, 18, 15, 26, 14, 27, 13, 17,
           24, 25, 21, 29, 12, 20, 28,  6,  7,  4, 19,  1,  3])

    '''
    n = np.size(x, 1)
    #initialize the array to return
    array = np.empty(n, dtype = np.int)
    
    #loop through rows
    for i in xrange(n):
        #loop through columns
        for j in xrange(n):
            #assign location in array if 1 is found
            if x[i,j] == 1:
                array[i] = j
                break
    return array
                            
def round_to_perm(x):
    '''
    Starting from a matrix close to a permutation matrix, round it to a
    permutation matrix.
    
    :param x: the "almost" permutation matrix
    :type x: a matrix
    :returns: a valid permutation matrix
    :rtype: a matix of floats
    
    Examples:
    
    >>> from pyQAP import *                                                   
    >>> import numpy as np
    >>> np.random.seed(123)
    >>> x = np.random.random((4,4))
    >>> x                          
    array([[ 0.69646919,  0.28613933,  0.22685145,  0.55131477],
           [ 0.71946897,  0.42310646,  0.9807642 ,  0.68482974],
           [ 0.4809319 ,  0.39211752,  0.34317802,  0.72904971],
           [ 0.43857224,  0.0596779 ,  0.39804426,  0.73799541]])
    >>> x = np.random.random((10,10))
    >>> round_to_perm(x)
    array([[ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])
    >>> x = np.random.random((13,13))
    >>> round_to_perm(x)
    array([[ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
           [ 0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.]])
    '''
    n = x.shape[0]
    p = np.zeros((n,n))
    for i in xrange(n):
        max_yet = -float('inf')
        for j in xrange(n):
            if x[i,j] > max_yet:
                max_yet = x[i,j]
                pos = j
        p[i,pos] = 1
    return p
    
def normalize(m):
    '''
    Return a normalized version of the matrix of size n. The new entries will
    be from 1 to (n^2-n)/2 with the exception of the diagonal which will remain
    filled with zeros.
    
    Examples:
    
    Unused and untested.
    '''
    n = len(m[0])
    tmp_list = []
    #iterate through and add top triangular entries to tmp_list
    for i in xrange(n):
        for j in xrange(i+1,n):
            tmp_list.append((m[i,j],i,j))
    
    #initiate new matrix to return
    ret = np.empty_like(m)
    
    #assign diagonal to zeros
    for i in xrange(n):
        ret[i,i] = 0.0
    
    #sort tmp_list
    get_key = lambda e: e[0]
    tmp_list.sort(key=get_key)
    
    #place integer versions of elements back into array
    current = 0
    for e in tmp_list:
        current +=1
        ret[e[1],e[2]] = current        #add main entry
        ret[e[2],e[1]] = current        #add symmetric entry
    
    return ret
